First-principles thermodynamics of transition metals: W, NiAl, PdTi 
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We apply the pseudopotential density functional perturbation theory approach along with the 
quasiharmonic approximation to calculate the thermal expansion of tungsten and two important 
metallic alloys, NiAl and PdTi. We derive the theory for anisotropic crystal structures and test the 
approximation that the anisotropic effects of thermal expansion are equivalent to negative pressure 
- this simplifies the calculation enormously for complex structures. Throughout, we find excellent 
agreement with experimental results. 
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I. INTRODUCTION 

First-principles calculations have established an ex- 
cellent record for describing the ground state proper- 
ties of materials at zero temperature over almost twenty 
years. More recently, the full description of the elec- 
tronic structure allows for accurate calculation of phonon 
frequencies|l| and the thermodynamic properties deriv- 
able from themfl- These methods have been applied to 
insulating systems and elemental metals. In this paper, 
we extend the methodology to investigate transition met- 
als and compounds that have a low symmetry or phonon 
anomalies. Previous authors have looked at elemental 
metals Al, Li, Na0, AgH, Cu@, and W@ and found 
that the quasiharmonic approximation holds most of the 
way to melting. We extend this work to alloys, which in- 
troduces some new features. As our prototype materials, 
we study W, NiAl and PdTi, which have many similari- 
ties but each of which has aspects of peculiar interest. 

NiAl is a high melting point alloy which finds appli- 
cations in aerospace and turbine blade design. Its elec- 
tronic structure is characterized by Kohn anomalies in 
the phonon spectrum 0, arising from a Fermi surface 
nesting effect. A huge amount of work has been done 
on the system using empirical potentials Ifl fiol ITU Il2| . 
which cannot describe the Kohn anomaly, and it is inter- 
esting to examine whether this neglect has serious struc- 
tural consequences. Crystallographically, NiAl has the 
cubic B2 (CsCl) crystal structure which is stable at all 
temperatures. This high symmetry simplifies the elec- 
tronic structure calculation enormously. 

PdTi provides a sharp contrast to NiAl. It is a high 
temperature shape memory alloy, isostructural with NiAl 
at high temperatures in the B2 phase, but undergoing a 
martensitic phase transition to the tetragonal B19 phase 
at 810 K and predicted to undergo a further transition to 
the monoclinic B19' at very low temperatures 01 ■ The 
dynamically stabilized high temperature phase and low 
symmetry low temperature phase each provide a chal- 
lenging test for theoretical prediction. 

Tungsten is another high melting point transition 
metal. In addition to studying its own thermal expan- 
sion, we use this to investigate the effects of anisotropic 
strain on the phonon contribution to the free energy. This 



enables us to make a simple test of some approximations 
required in the calculation of PdTi. 

The difference between the materials being studied is 
highlighted by the phonon spectra for the B2 phase: NiAl 
and PdTi both adopt the B2 structure at high temper- 
ature while tungsten is a refractory metal and has a bcc 
structure at all temperature. The bcc structure can be 
regarded as a special case of B2 where both atoms are the 
same. NiAl has phonon anomalies due to Fermi surface 
nesting, while B2 structure PdTi is dynamically stabi- 
lized: it has negative frequency phonons (i.e. mechan- 
ically unstable) at low temperatures. Such a situation 
leads to temperature dependent modes, which cannot 
be treated within the quasiharmonic approximation and 
would need to be dealt with separately either in real[l4| 
or reciprocal ^ space. 

The quasiharmonic approximation^^ to the free en- 
ergy assumes all phonons can be treated as simple har- 
monic oscillators with a frequency dependent on the vol- 
ume of the material. The static lattice energy (T = 0) is 
evaluated at a range of volumes, and for each of these the 
phonon dispersion relation is calculated. The quasihar- 
monic contribution to the free energy at finite temper- 
ature is then evaluated from Bose-Einstein statistics of 
the phonons at a given volume. From this the pressure 
can be evaluated for all conditions of T and V, giving 
an equation of state. Alternately in the zero pressure 
case, the relationship between T and V (i.e. the ther- 
mal expansion) can be determined. Since larger volumes 
typically lead to lower phonon frequencies, the entropy is 
typically lowered by expansion. Thus larger volume (en- 
tropy) is favoured at higher temperatures, and thermal 
expansion occurs. 

There are, of course, other contributions to thermal 
expansion in ordered alloys (antisite defect creation, va- 
cancy creation, coupling of phonons to lattice parame- 
ters, anharmonic effects due to finite atomic displace- 
ments) and it is interesting to test how much of the ob- 
served expansion can be attributed to quasiharmonic ef- 
fects. 
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II. THEORY 

Phonon dispersions can be evaluated using density- 
functional perturbation theory (DFPT) [nj and finite 
displacement methods[l8|. Here we use the former. Pre- 
vious work on isotropic equations of state for cubic ma- 
terials covers much of the theory, here we investigate the 
practical requirements for lower symmetry materials, in 
particular what approximations can be made to made 
calculations tractable. 

Virtually all implementations of density-functional 
theory (DFT) based ab initio codes take as input the unit 
cell vectors a, b and c and atomic positions. We define 
a matrix B comprised of the three unit cell vectors. 

B = (a, b, c) 

The volume of the unit cell is 

V = \B\ = (a x b) .c. 

For a cubic material, the quasiharmonic approxima- 
tion uses the phonon frequencies calculated at different 
volume to obtain the free energy as a function of two 
variables, 

F(V,T) = E(V) - k B T In Z(V) (1) 

where E(V) is the cold curve energy at the volume 
V (i.e. the energy of a system with the atoms placed on 
their lattice sites), k B is the Boltzmann's constant and Z 
is the vibrational partition function. We assume that for 
all temperatures of interest the contribution of electronic 
excitations is negligible - generally true for equation of 
states [19|, but not for transport properties such as heat 
conduction. 

From F(V, T) one can obtain the equation of state of 
the single phase of the system: 

and the volume thermal expansion: 

aV= v{sT) p (3) 

These can be evaluated for a given structure regardless 
of whether an alternate phase of lower Gibbs free energy 
exists. 

For non-cubic materials, the free energy is a function 
of all independent cell parameters, and the linear ther- 
mal expansion is a tensor quantity. Calculating the free 
energy across such a multidimensional space is impracti- 
cal, so we examine a further approximation: that the ef- 
fect of increasing temperature is equivalent to a negative 



pressure[20| At each pressure, the internal coordinates 
must be re-relaxed to their equilibrium position within 
the unit cell, since the phonon spectrum is only well de- 
fined by an expansion about such an equilibrium. 
Thus, we proceed as follows: 

1/ Evaluate the equilibrium structure at P = 0, by 
minimizing f7(a, b, c, {u;}) = Uo where Ui are the inter- 
nal atomic coordinates. 

2/ Evaluate the force constants and hence the phonon 
dispersion relation for the relaxed structure, and the par- 
tition function Z(V) for these oscillators. V — \B\ = 
(a x b) . c. 

3 / Evaluate Gibbs free energy at a range of tempera- 
tures F(V, T) = U - fcTln Z + PV. 

4/ Evaluate the equilibrium structure at P, by mini- 
mizing enthalpy U (a, b, c, {ui}) + PV = H (V) where -Uj 
are the internal atomic coordinates. 

5 / Evaluate the phonon dispersion relation for this re- 
laxed structure, and the partition function Z for these 
oscillators. 

6/ Evaluate Gibbs free energy at fixed V for a range 
of temperatures F(V, T) = U - kT In Z + PV. 

7/ Repeat 4, 5, 6 for a range of positive and negative 
pressures. 

Once a coarse grid of 1^-points has been generated (for 
thermal expansion, as few as 3 different pressures suffices) 
a finer grid can be generated without further expensive 
ab initio calculation by interpolating force constants ^ij- 
A dense grid of T-points can easily be generated at a 
given V. 

In terms of speed and accuracy there is little to choose 
between DFPT or finite displacements: the former in- 
volves calculation of second derivatives and the later re- 
quires large supercells. We note that for finite displace- 
ments, once the eigenvectors are calculated at P = 0, a 
single calculation of restoring forces with all atoms dis- 
placed suffices at other volumes @- 

The allowed occupation states of each phonon mode of 
frequency u> are e n (u>) = (n + 1/2) fko where n is the the 
number of phonons populating that mode. The canonical 
partition function is therefore 



= e- pe "^ = l/(2sinh(^/2fc B r)). (4) 

71 = 

The total partition function is then: 

z = ru( w ) (5) 

For calculating thermal expansion under constant 
(zero) pressure boundary condition, the essential quan- 
tity is the free energy. The mean vibrational free energy: 



-k B T In Z(V) = -k B T 9 {cj,V) ln(£(w, T)) du, (6) 
Jo 



3 



introducing the phonon density of states g(oj, V) to 
make the dependence of density of states on volume ex- 
plicit. 

Phonon frequencies are evaluated at a dense set of q- 
points in the first Brillouin zone using the force constants 
matrix obtained by fourier transform of the "exact" 
DFPT frequencies. The density of phonon states is cal- 
culated by integration using the tetrahedron method [2lj. 
This provides the Z(V) required in step 2/ above. 

The structure obtained at a particular temperature can 
be found by minimizing the free energy with respect to 
the lattice vectors and internal coordinates of the atoms 



dF 
dBi4 



= 0; 



dF 



= 



(7) 



These derivatives can be split into cold curve 
(E(B, Uk)) and phonon (F p hon) contributions. 

To make calculations for complex structures tractable, 
we make the assumption that the phonon part of the 
free energy is a function of volume only, and is indepen- 
dent of Uk- We examine this approximation in the case 
of anisotropic deformation of tungsten and find it to be 
good. Intuitively, this can be understood if the force 
constants vary linearly with separation: we might expect 
that volume conserving shears and changes in Uk will 
stiffen some force constants and weaken others, raising 
some frequencies and lowering others to give an over- 
all cancellation in the small strain limit. By contrast, 
volume increases weaken all force constants and give a 
systematic decrease in mean frequency. Hence the con- 
ditions above become: 



dE 

dB~ W 



dF, 



phon I 



dv |Uk 

dE 



(8) 

(9) 
(10) 



III. RESULTS 



1. Anisotropic strain - Tungsten 



We calculated phonon density of states (DOS) at 15 
different structures corresponding to tetragonal strains 
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dispersion of bec tungsten at equilibrium lattice parame- 
ter along with experimental results is shown in Fig. ^ 
while the phonon dispersions under tetragonal strains 
and a larger volume are shown in Fig.|2Ia) and Fig.^b), 
respectively. The calculated phonon frequencies are in 
good agreement with experiment (Fig. ^) and with a re- 
cent calculation §|. Fig.^a) and (b) show that compared 
to volumetric strain, volume-conserving shear strains 
have negligible contribution to phonon free energy. 

The free energy was evaluated at each strain state. In 
Fig. |3 we show the thermal expansion calculated using 
the three isotropic strains. The linear thermal expansion 
adopts the form from Ref. The agreement between 
theoretical and experimental results is very good. 

Fig. 01 shows the phonon free energy with different 
strain. As expected, the phonon contribution to the 
free energy is only weakly affected by volume-conserving 
anisotropic strain. Thus it is a good approximation to 
neglect this term: a given volume-conserving strain con- 
tributes about 5% as much as equivalent volume strain. 



2. Anomalous phonons - NiAl 

The phonon spectrum of B2 NiAl is characterized by 
Kohn anomalies. As a consequence, to describe the 
phonon dispersion properly requires a thorough DFPT 
q-point sampling or (equivalcntly) large supercclls in di- 
rect space methods. In Fig. [5] we show the convergence 
of the result to experimental observation with (/-point 
sampling density. 

Despite the anomalies, we find that the thermody- 
namic properties are almost independent of g-point sam- 
pling aside from a rigid shift of the free energy. This 
does not affect the thermal expansion, but stabilizes the 
structure slightly against alternate structures because of 
the increased entropy of the low frequency phonons. 

The entire calculated phonon dispersion relations at 
two different volumes are shown in Fig.EJa). All phonon 
frequencies increase under compression. Hence the mode 
Griineisen parameters, 7j(q), shown in Fig. EJb), and 
defined as 



We use tungsten as a benchmark to investigate the ef- 
fects of anisotropic strain on the phonon contribution to 
the free energy. Tungsten has a bec crystal structure. 
The calculations were done with the generalized gradi- 
ent approximation (GGA) and the equilibrium lattice pa- 
rameter was computed to be 2.905 A, comparable to the 
experimental result of 2.887 A. Phonon calculations were 
done at three different volumes: Vb, 1.01 3 Vo and 1.02 3 Vb. 
At each volume, phonon dispersion was done with the bec 
structure and tetragonal strained structures with four dif- 
ferent c/a ratios, i.e. 0.96, 0.98, 1.02 and 1.04. 



are positive throughout the whole BZ. There is no 
anomalous thermal expansion. Integration of free ener- 
gies derived from phonon spectra gives the thermal ex- 
pansion shown in Fig. [7| The excellent agreement with 
experimental data suggests that both DFPT phonons 
and the quasiharmonic approximation are valid in this 
case. 
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3. Complex structure - PdTi 

For PdTi we calculate phonon spectra for two phases, 
the low temperature B19 1 and the ambient tempera- 
ture B19 phase. In both cases, the equilibrium cold 
curve structure is determined for a particular pressure 
by relaxing internal parameters and lattice constants 
simultaneously It is then assumed that the effect 
of thermal expansion on lattice parameter ratios and 
internal parameters is equivalent to that of (negative) 
pressure, such that the phonon spectrum need only be 
evaluated once at each volume, and the free energy can 
be interpolated as a function of volume only, with the 
free parameter relaxation done on the cold-curve struc- 
ture. The results for anisotropic deformation in tungsten 
(Fig. suggest that even if the effects of pressure and 
thermal expansion on free parameters are not equivalent, 
the consequent change in thermal free energy due to non- 
volumetric strains will be small. 

The B19' phase has stable phonons throughout the 
entire BZ (Fig.^a)) but the B19 phase (Fig.^b)) has 
some imaginary frequency modes: we assume these con- 
tribute to the free energy like free particles [l9j. Assuming 
the phonons in the B19' phase to remain harmonic up to 
and above the phase transition, this enables us to cal- 
culate the relative free energies of the two phases, and 
determine for each temperature not only the lattice pa- 
rameters but also which is the stable phase. This gives 
rise to the thermal expansion shown in Fig. [5] with a dis- 
continuity marking the phase transition. 

The calculated phase transition temperature for the 
B19-B19' phase, on the basis of equal quasiharmonic 
free energies, is 140K. This can be compared with an- 
other estimate based on treating the 519 phase as a 
barrier between equivalent B19' variants, which barrier 
(about 0.0007eV/atom) can continually be crossed when 
the temperature reaches AE/ks = 9K in the higher tem- 
perature phase. The discrepancy between these arises be- 
cause one concentrates on dynamics of a particular soft 
mode, while the other considers the static average over all 
modes: experimentally, the lower of the two is expected. 



mation and DFPT finding that agreement with experi- 
ment is excellent. This shows that the harmonic phonon 
free energy is by far the dominant effect in thermal ex- 
pansion, and that other effects such as thermal defects, 
electron free energy, coupling of phonons to lattice pa- 
rameters, and anharmonic effects due to finite atomic 
displacements can safely be neglected. 



These three materials are chosen to represent repre- 
sent an increasingly challenging test to this methodol- 
ogy. The accurate thermal expansion in tungsten was 
expected on the basis of previous work [HQ, and allowed 
for the demonstration that anisotropic strain has only 
a second order effect on the free energy. Similarly the 
results on NiAl show that the method can be applied to 
alloys just as effectively as with elements and showed that 
the presence of phonon anomalies also has little effect on 
thermal expansion. 



Without these two results, the PdTi calculation would 
have been intractable, but armed with the knowledge 
that the anisotropic temperature effect can be treated 
as a pressure effect on the cold curve, allowing the vi- 
brational partition function to be treated as a func- 
tion of volume only, calculation and minimizations of 
the full F(By,Uk,T) becomes tractable. In addition, 
implementing the previously-described treatment of soft 
phonons 0, [h| allow us to predict the thermal expan- 
sion of both B19' and the dynamically stabilized B19 
PdTi crystal, in addition to estimating the phase transi- 
tion temperature for B19-B19'. To our knowledge, these 
quantities have yet to be measured experimentally and 
such measurement will provide a sensitive test to our 
methods. 



In sum, this paper represents a significant theoretical 
advance in the type of materials whose thermal expansion 
can be calculated from ab initio simulation. 



IV. DISCUSSION 

We have calculated the thermal expansion properties 
of W, NiAl and PdTi using the quasiharmonic approxi- 
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FIG. 1: Phonon spectra calculated for bcc W at ao=2.905A, 
using the PWSCF and PHONON code with DFPTlH 
and a norm conserving pseudopotential. The electronic wave 
functions were represented in a plane-wave basis set with a 
kinetic energy cutoff of 32 Ry. The Brillouin zone (BZ) inte- 
grations were carried out by the Hermite-Gaussian smearing 
technique [25l using a 16 x 16 x 16 fc-point mesh. Phonon calcu- 
lations are exact on meshes which are commensurate with the 
g-point mesh: a 8 x 8 x 8 mesh for tungsten. Other phonons 
come from dynamical matrices calculated using force con- 
stants determined by Fourier transform of original data 24] . 
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FIG. 2: [001] and [100] branches of the tungsten dispersion 
relation for volume-conserving strains of (a) 4% on (001) un- 
strained and (b) bulk strains of 1% and 2%. The key obser- 
vation is that for volume conserving shear some frequencies 
increase while others decrease, meanwhile for volume expan- 
sion all frequencies decrease. Thus in the former case some 
compensation occurs to keep the phonon free energy constant, 
while in the latter case all modes contribute to increase free 
energy. Calculation details are as for Fig. 
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FIG. 3: Calculated linear thermal expansion for tungsten, 
compared with experimental data from Refl26t ao represents 
the room temperature volume. 
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FIG. 4: Graphs of tungsten phonon free energy against tem- 
perature for 15 different strain conditions, plotted relative to 
G(V) for V at the cold curve minimum. The effect of volume 
strain can be seen to be about ten times greater than that of 
shear strain. 
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FIG. 5: Phonon dispersion relations along the acoustic (110) 
branches in NiAl, calculated with increasing (/-point density. 
Lines are calculated data, taken from force constants deduced 
from Fourier transforming the ui(q) data. Labels indicate the 
number of q-points used in the DFPT calculation. Symbols in- 
dicate experimental data taken from Ref. for a composition 
of Ni 50 Al 5 o Calculations use the PWSCF and PHONON||3] 
code and DFPT |2^I using an ultrasoft pseudopotentials for 
Ni, |2Sj. and a norm conserving pseudopotential for Al. The 
electronic wave functions were represented in a plane-wave 
basis set with a kinetic energy cutoff of 30 Ry. The BZ 
integrations|2fil| used 12 x 12 x 12 mesh for NiAl, giving exact 
phonon on a 6 x 6 x 6 mesh. 
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FIG. 6: (a) Phonon spectra calculated for B2 NiAl, solid 
lines are at a lattice parameter of ao=2.906A(GGA theoreti- 
cal lattice constant) the dotted lines at a lattice of 1.02oo. All 
phonons have higher frequencies under compression, (b) Cal- 
culated mode Griineisen parameters of NiAl along symmetry 
lines of simple cubic BZ at equilibrium volume. 
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FIG. 7: Calculated thermal expansion for NiAl, compared 
with experimental data from Ref. I29I The zero of lattice pa- 
rameter ao is chosen to be 2.905Ain the theory and 2.886Ain 
the experimental case, which correspond to 273K 
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FIG. 8: Phonon spectra calculated for (a) B19' PdTi 
(b) B19 PdTi. B19 has many soft phonon modes which 
characterize the dynamical instability (grey region, these 
frequencies are actually imaginary). Calculations ultrasoft 
pseudopotentials^3 12 x 8 x 8 mesh and 30Ry cutoff. 
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FIG. 9: Calculated percentage thermal expansion for PdTi, 
including the B19-B19' phase transition (expanded scale in 
lower insert). The dotted lines shows the calculated thermal 
expansion for the thermodynamically unstable phase (B19 at 
low T, B19' at high T), assuming that all modes remain har- 
monic. The reference volume is the T=0 volume for B19': 
30.64A 3 which is indistinguishable from the equilibrium for 
B19: 30.64A 3 and about 0.34% larger than the cold curve 
minimum[13j which excludes zero-point vibrations. The up- 
per insert shows the free energy difference between the two 
phases as a function of temperature, with the phase transi- 
tion at 140K. This is very much larger than the estimate taken 
from the Landau barrier height [l3| 



